# How many cases and controls are there?
sum(CaseControl$case, na.rm=T)
sum(CaseControl$case==0, na.rm=T)
coordinates(CaseControl) = ~long + lat
class(CaseControl)  # class is now SpatialPointsDataFrame
CaseControl@data # View attribute data
# Visualise cases and controls
plot(StudyArea)
plot(CaseControl[CaseControl@data$case==0,], col="blue", pch=19, add=TRUE)
plot(CaseControl[CaseControl@data$case==1,], col="red", pch=19, add=TRUE)
# Take cases only
cases<-CaseControl[CaseControl@data$case==1,]
controls<-CaseControl[CaseControl@data$case==0,]
coordinates(cases)
CasesPPP<-as(cases, "ppp")
# Ripley's K function applied only to cases
K<-Kest(CasesPPP)
plot(K)
E<-envelope(CasesPPP, Kest, nsim=999)
plot(E)
KX <- Kest(CaseControlPPP[CaseControlPPP$marks==1])
plot(KX, sqrt(iso/pi) ~ r)
coordinates(CaseControl) = ~long + lat
class(CaseControl)  # class is now SpatialPointsDataFrame
CaseControl@data # View attribute data
CaseControlPPP<-as(CaseControl, "ppp")
KX <- Kest(CaseControlPPP[CaseControlPPP$marks==1])
CaseControlPPP<-ppp(CaseControl$long, CaseControl$lat, range(CaseControl$long), range(CaseControl$lat), marks = as.factor(CaseControl$IndexCaseHousehold))
KX <- Kest(CaseControlPPP[CaseControlPPP$marks==1])
CaseControlPPP<-ppp(CaseControl$long, CaseControl$lat, range(CaseControl$long), range(CaseControl$lat), marks = as.factor(CaseControl$case))
KX <- Kest(CaseControlPPP[CaseControlPPP$marks==1])
plot(KX, sqrt(iso/pi) ~ r)
KY <- Kest(CaseControlPPP[CaseControlPPP$marks==0])
plot(KY, sqrt(iso/pi) ~ r)
Kdiff <- eval.fv(KX - KY)
plot(Kdiff)
kenv.label(Kdiff) ## says this function exists in Bivand but not found
library(smacpod) # statistical methods for the analysis of case-control point data
length(levels(CaseControlPPP$marks))
levels(CaseControlPPP$marks)
kdest = kdest(CaseControlPPP, case = 2,nsim=999) # more simulations = more significance
kdest
plot(kdest)
kdplus.test(kdest)
?kdest
kdest = kdest(CaseControlPPP, case = 2,nsim=999, level=0.05) # more simulations = more significance
kdest = kdest(CaseControlPPP, case = 2,nsim=999, level=0.95) # more simulations = more significance
plot(kdest)
kdplus.test(kdest)
?plot
??plot
? kdplus.test
# Clean workspace
rm(list = ls())
# Load paths to file directories (adapt as needed)
root<-"C://"
path<-"Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//JLS - Point processes//"
# Attach libraries for visualisation
library(rgdal)
library(geoR)
# Attach libraries for point processes
library(raster)
# Attach library to convert between data types used by different packages (ppp and spatial points)
library(sp)
library(maptools) # Must have rgeos installed as well
###########################################################################################################
path<-"Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//Lecture 5 Analysis of spatial clustering//"
setwd(paste(root, path, "Lab 1", sep=""))
load("GambiaID.RData")
Gambia <- readRDS("GambiaID.RData")
GMB_Adm2<-getData('GADM', country='GMB', level=2)
proj4string(GMB_Adm2) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ')
# Calculate prevalence from individual malaria data by cluster ID
GambiaPrev<-as.data.frame(cbind(unique(Gambia$x), unique(Gambia$y), tapply(Gambia$pos, Gambia$ID, mean)))
colnames(GambiaPrev)<-c("x", "y", "Prev")
map <- fortify(GMB_Adm2, region="NAME_2")                                                  #"ggplot2" package
ggplot() +
geom_path(data = map, aes(x = long, y = lat, group = group)) +
geom_point(data = GambiaPrev, aes(x = x, y = y, size = Prev), color = "red")
library(ggplot2)
map <- fortify(GMB_Adm2, region="NAME_2")                                                  #"ggplot2" package
ggplot() +
geom_path(data = map, aes(x = long, y = lat, group = group)) +
geom_point(data = GambiaPrev, aes(x = x, y = y, size = Prev), color = "red")
coordinates(GambiaPrev)<-~x+y # convert to SPDF                                            # "sp" package
# adding Coordinate Referent Sys.
proj4string(GambiaPrev) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ')
gambia.dists <- as.matrix(dist(cbind(GambiaPrev$x, GambiaPrev$y)))
# Take the inverse of the matrix values so that closer values have a larger weight and vs vs
gambia.dists.inv <- 1/gambia.dists
diag(gambia.dists.inv) <- 0 # replace the diagonal values with zero
?Moran.I
??Moran.I
Moran.I(GambiaPrev$Prev, gambia.dists.inv) # Global Moran's I test                           # "ape" package
library(ape) # Analyses of Phylogenetics and Evolution
Moran.I(GambiaPrev$Prev, gambia.dists.inv) # Global Moran's I test                           # "ape" package
maxDist<-max(dist(cbind(GambiaPrev$x, GambiaPrev$y)))
maxDist
xy=cbind(GambiaPrev$x, GambiaPrev$y)
pgi.cor <- correlog(coords=xy, z=GambiaPrev$Prev, method="Moran", nbclass=10)
library(pgirmess) # Data Analysis in Ecology
pgi.cor <- correlog(coords=xy, z=GambiaPrev$Prev, method="Moran", nbclass=10)               # "pgirmess" package
# coordys = xy corredinates, z= vector of values at each location and nbclass = the number of bins
plot(pgi.cor) # statistically significant values (p<0.05) are plotted in red
Moran.I(GambiaPrev$Prev, gambia.dists.inv)                                                   # "ape" package
maxDist<-max(dist(cbind(GambiaPrev$x, GambiaPrev$y)))
maxDist
xy=cbind(GambiaPrev$x, GambiaPrev$y)
pgi.cor <- correlog(coords=xy, z=GambiaPrev$Prev, method="Moran", nbclass=10)               # "pgirmess" package
coords<-coordinates(xy)
class(coords)
IDs<-row.names(as.data.frame(coords))
# In this approach, we chose a distance d such that pairs of points with distances less than
# d are neighbors and those further apart are not.
Neigh_nb<-knn2nb(knearneigh(coords, k=1, longlat = TRUE), row.names=IDs)
#assigns at least on neighbor to each and calculates the distances between
dsts<-unlist(nbdists(Neigh_nb,coords))
# returns the distance between nearest neighbors for each point
summary(dsts)
??knn2nb
library(spdep) # Spatial Dependence: Weighting Schemes, Statistics and Models
Neigh_nb<-knn2nb(knearneigh(coords, k=1, longlat = TRUE), row.names=IDs)                 # "spdep" package
#assigns at least on neighbor to each and calculates the distances between
dsts<-unlist(nbdists(Neigh_nb,coords))
# returns the distance between nearest neighbors for each point
summary(dsts)
max_1nn<-max(dsts)
max_1nn # maximum distance to provide at least one neighbor to each point
Neigh_kd1<-dnearneigh(coords,d1=0, d2=max_1nn, row.names=IDs)   # neighbors within maximum distance
Neigh_kd2<-dnearneigh(coords,d1=0, d2=2*max_1nn, row.names=IDs) # neighbors within 2X maximum distance
nb_1<-list(d1=Neigh_kd1, d2=Neigh_kd2)
nb_1<-list(d1=Neigh_kd1, d2=Neigh_kd2)
sapply(nb_1, function(x) is.symmetric.nb(x, verbose=F, force=T))
# Checks for symmetry (i.e. if i is a neighbor of j, then j is a neighbor of i). Does not always hold for k-nearest neighbours
sapply(nb_1, function(x) n.comp.nb(x)$nc)
par(mfrow=c(1,1))
plot(xy, pch=16)
plot(Neigh_kd1, coords, col="green",add=T)
plot(xy, pch=16)
plot(Neigh_kd2, coords,col="green", add=T)
weights<-nb2listw(Neigh_kd1, style="B")   # simplest binary weights, using minimum distance for one neighbor
weights
sum(weights)
weights<-nb2listw(Neigh_kd1, style="B")   # simplest binary weights, using minimum distance for one neighbor
weights
##STEP 3
moran.test(GambiaPrev$Prev, listw=weights)  #using row standardised weights
set.seed(1234)
bperm<-moran.mc(GambiaPrev$Prev, listw=weights,nsim=999)
bperm
#statistic = 0.582, observed rank = 1000, p-value = 0.001
# Plot simulated test statistics
par(mfrow=c(1,1))
hist(bperm$res, freq=T, breaks=20, xlab="Simulated Moran's I")
abline(v=0.58284, col="red")
I <-localmoran(GambiaPrev$Prev, weights)
Coef<-printCoefmat(data.frame(I[IDs,], row.names=row.names(coords),
check.names=FALSE))
LocalI<-as.data.frame(cbind(coords,Coef))
LocalI$C_mI <- LocalI[,3] - mean(LocalI[,3])
nci<-moran.plot(GambiaPrev$Prev, listw=weights,
xlab="Standardised Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
GambiaPrev$Prev
GambiaPrev$StdPrev<-(GambiaPrev$Prev-mean(GambiaPrev$Prev))/sd(GambiaPrev$Prev)
nci<-moran.plot(GambiaPrev$Prev, listw=weights,
xlab="Standardised Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
nci2<-moran.plot(GambiaPrev$StdPrev, listw=weights,
xlab="Standardised Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
nci<-moran.plot(GambiaPrev$Prev, listw=weights,
xlab="Standardised Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
nci2<-moran.plot(GambiaPrev$StdPrev, listw=weights,
xlab="Standardised Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
weights <- lag.listw(weights,GambiaPrev$StdPrev)
weights2 <- lag.listw(weights,GambiaPrev$StdPrev)
weights<-nb2listw(Neigh_kd1, style="B")   # simplest binary weights, using minimum distance for one neighbor
weights
weights2 <- lag.listw(weights,GambiaPrev$StdPrev)
nci2<-moran.plot(GambiaPrev$StdPrev, listw=weights2,
xlab="Standardised Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
morlm <- lm(weights2 ~ GambiaPrev$StdPrev)
aa <- morlm$coefficients[1]
mori <- morlm$coefficients[2]
par(pty="s") # make sure it’s square
par(pty="s") # make sure it’s square
plot(GambiaPrev$StdPrev,weights,xlab="Standardised Prevalence",ylab="Spatial Lag of Prevalence")
plot(GambiaPrev$StdPrev,weights2,xlab="Standardised Prevalence",ylab="Spatial Lag of Prevalence")
abline(aa,mori)
abline(h=0,lty=2)
abline(v=0,lty=2)
title(paste("Moran Scatterplot I= ",format(round(mori,4))))
nci<-moran.plot(GambiaPrev$Prev, listw=weights,
xlab="Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
Coef
I
moran.test(GambiaPrev$Prev, listw=weights)  #using row standardised weights
require(GeoXp)
#map local outliers
infl<-apply(nci$is.inf,1,any)
sum(infl==T)#7 true (11% - more than would expect by chance)
x<-GambiaPrev$Prev
lhx<-cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=T)
wx<-lag(W.list,GambiaPrev$Prev)
lhwx<-cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=T)
lhlh<-interaction(lhx,lhwx,infl,drop=T)
cols<-rep(1, length(lhlh))
cols[lhlh=="L.L.TRUE"]<-2
cols[lhlh=="H.L.TRUE"]<-3
cols[lhlh=="L.H.TRUE"]<-4
cols[lhlh=="H.H.TRUE"]<-5
sum(cols>1) ##7 outliers
# Map points that are local outliers in the plot
infl<-apply(nci$is.inf,1,any)
sum(infl==T)#7 true (11% - more than would expect by chance)
nci
infl
sum(infl==T)#7 true (11% - more than would expect by chance)
x<-GambiaPrev$Prev
lhx<-cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=T)
wx<-lag(weights,GambiaPrev$Prev)
lhwx<-cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=T)
lhlh<-interaction(lhx,lhwx,infl,drop=T)
cols<-rep(1, length(lhlh))
cols[lhlh=="L.L.TRUE"]<-2
cols[lhlh=="H.L.TRUE"]<-3
cols[lhlh=="L.H.TRUE"]<-4
cols[lhlh=="H.H.TRUE"]<-5
sum(cols>1) ##7 outliers
plot(GMB_Adm2,border="darkgrey")
points(xy, col=c("lightgrey", "cyan", "coral","coral4","cyan4")[cols], pch=16, cex=0.6, add=T)
legend("topright", legend=c("None", "LL","HL", "LH", "HH"),
fill=c("lightgrey", "cyan", "coral","coral4","cyan4"), bty="n", cex=0.8,
y.intersp=0.8)
nci<-moran.plot(GambiaPrev$Prev, listw=weights,
xlab="Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")
rm(list = ls())
# Load paths to file directories (adapt as needed)
root<-"C://"
path<-"Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//JLS - Point processes//"
# Attach libraries for visualisation
library(rgdal)
library(geoR)
library(gcmr)
# Attach libraries for point processes
library(sp)
library(spatstat)
library(smacpod) # Spatial scanning statistic
# Attach library to convert between data types used by different packages (ppp and spatial points)
library(splancs)
library(maptools) # Must have rgeos installed as well
library(plotrix) #  Plotting spatial scanning
###########################################################################################################
#  Point processes
###########################################################################################################
# 1. Simulating complete spatial randomness (CSR): different realisations of a
# random spatial process will appear different but have same underlying properties
# Waller, L.A. and Gotway, C.A. (2004) Applied Spatial Statistics
#   for Public Health Data.  John Wiley and Sons: New York.
par(mfrow=c(2,3),pty="s")
for (i in 1:6) {
# x and y coordinates are uniform random variates (see text for discussion)
x <- runif(30) #generates 30 random deviates from the uniform distribution
y <- runif(30)
plot(x,y,xlim=c(0,1),ylim=c(0,1),xlab="u",ylab="v",cex.lab=1.5,cex.axis=1.1)
}
###########################################################################################################
#  Part II:  Global cluster detection for event data
###########################################################################################################
root<-"C://"
path<-"Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//Lecture 5 Analysis of spatial clustering//"
setwd(paste(root, path, "Lab 2", sep=""))
# Read in case-control data
CaseControl<-read.csv("CaseControl_v4.csv", header=T, stringsAsFactors=FALSE)
CaseControl<-read.csv("CaseControl.csv", header=T, stringsAsFactors=FALSE)
project2<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 " #
StudyArea<-readShapePoly(paste(root, path, map.dir, "StudySiteDistricts.shp", sep=""),IDvar=NULL, proj4string=CRS(project2))
StudyArea<-readShapePoly("StudySiteDistricts.shp",IDvar=NULL, proj4string=CRS(project2))
sum(CaseControl$case, na.rm=T)
sum(CaseControl$case==0, na.rm=T)
# Promote the case data to a PPP data type
CasesPPP<-as(cases, "ppp")
cases<-CaseControl[CaseControl@data$case==1,]
CaseControl<-read.csv("CaseControl.csv", header=T, stringsAsFactors=FALSE)
project2<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 " #
StudyArea<-readShapePoly("StudySiteDistricts.shp",IDvar=NULL, proj4string=CRS(project2))
###           ADD CODE    #2                     ##
# Plot the data using techniques you have already learned
# Q1: Visually assess the presence of clustering across the study area.
# Promote data frame to spatial
coordinates(CaseControl) = ~long + lat
class(CaseControl)  # class is now SpatialPointsDataFrame
CaseControl@data # View attribute data
# Take cases only
cases<-CaseControl[CaseControl@data$case==1,]
controls<-CaseControl[CaseControl@data$case==0,]
coordinates(cases)
# back to data
#CaseControl<-as.data.frame(CaseControl)
###########################################################################################################
#  Ripley's K function
###########################################################################################################
# Promote the case data to a PPP data type
CasesPPP<-as(cases, "ppp")
K<-Kest(CasesPPP)
plot(K)
E<-envelope(CasesPPP, Kest, nsim=999)
plot(E)
par(mfrow=c(1,1))
plot(K)
E<-envelope(CasesPPP, Kest, nsim=999)
plot(E)
# Create a marked point process
CaseControlPPP<-ppp(CaseControl$long, CaseControl$lat, range(CaseControl$long), range(CaseControl$lat), marks = as.factor(CaseControl$case))
# Calculate the K-function for cases
KX <- Kest(CaseControlPPP[CaseControlPPP$marks==1])
plot(KX, sqrt(iso/pi) ~ r)
# Calculate the K-function for controls
KY <- Kest(CaseControlPPP[CaseControlPPP$marks==0])
plot(KY, sqrt(iso/pi) ~ r)
Kdiff <- eval.fv(KX - KY)
plot(Kdiff)
kenv.label(Kdiff) ## says this function exists in Bivand but not found
?plot
plot(Kdiff, legendpos="float")
KE <- envelope(CaseControlPPP, Kdiff, nsim=19)
plot(K)
length(levels(CaseControlPPP$marks))
levels(CaseControlPPP$marks)
kdest = kdest(CaseControlPPP, case = 2,nsim=999, level=0.95)                                        #"smacpod" package
kdest
plot(kdest)
kdplus.test(kdest)
sG<-Sobj_SpatialGrid(StudyArea)$SG
gt<-slot(sG,"grid")
# compute the risk ratio by estimating the intensity of cases and controls and taking the ratio
grid<-slot(slot(slot(StudyArea, "polygons")[[1]], "Polygons")[[1]], "coords")
mean(dist(coordinates(cases)))
grid
dist(unique(grid))
grid<-slot(slot(slot(StudyArea, "polygons")[[1]], "Polygons")[[1]], "coords")
mean(dist(coordinates(cases)))
s<-seq(0,0.12, by=0.0001) #A vector of distances at which to calculate the K function
khcases<-khat(coordinates(cases), grid, s)
khcontrols<-khat(coordinates(controls), grid, s)
# Covariance matrix (?)
khcov<-khvmat(coordinates(cases), coordinates(controls), grid, s)
khcov<-khvmat(coordinates(cases), coordinates(controls), grid, s)
# Test statistic
T0<-sum(((khcases-khcontrols))/sqrt(diag(khcov+0.00000000001)))
niter<-999  # Number of iterations to simulate
T<-rep(NA,niter) # Empty vector in which to store simulated values of T
khcasesrel<-matrix(NA, nrow=length(s), ncol=niter)
khcontrolsrel<-matrix(NA, nrow=length(s), ncol=niter)
# Simulate values
DIFF<-NULL
for(i in 1:niter){
idxrel<-sample(CaseControl$IndexCaseHousehold)==0
controlsrel<-coordinates(CaseControl[!idxrel,])
casesrel<-coordinates(CaseControl[idxrel,])
khcasesrel[,i]<-khat(casesrel, grid, s)
khcontrolsrel[,i]<-khat(controlsrel, grid, s)
DIFF<-cbind(DIFF,khdiff)
khdiff<-khcasesrel[,i]-khcontrolsrel[,i]
length(khdiff); length(sqrt(diag(khcov)))
T[i]<-sum(khdiff/sqrt(diag(khcov+0.00000000001)))  # check this line, did not work as entered
}
# Calculate the P-value (proportion of times the simulated value is more extreme than observed)
niter<-999  # Number of iterations to simulate
T<-rep(NA,niter) # Empty vector in which to store simulated values of T
khcasesrel<-matrix(NA, nrow=length(s), ncol=niter)
khcontrolsrel<-matrix(NA, nrow=length(s), ncol=niter)
# Simulate values
DIFF<-NULL
for(i in 1:niter){
idxrel<-sample(CaseControl$IndexCaseHousehold)==0
casesrel<-coordinates(CaseControl[idxrel,])
controlsrel<-coordinates(CaseControl[!idxrel,])
khcasesrel[,i]<-khat(casesrel, grid, s)
khcontrolsrel[,i]<-khat(controlsrel, grid, s)
khdiff<-khcasesrel[,i]-khcontrolsrel[,i]
DIFF<-cbind(DIFF,khdiff)
length(khdiff); length(sqrt(diag(khcov)))
T[i]<-sum(khdiff/sqrt(diag(khcov+0.00000000001)))  # check this line, did not work as entered
}
pvalue<-(sum(T>T0)+1)/(niter+1)
pvalue
# Estimate 95% CI by taking the maximum and minimum simulated value for each point estimated
UppEnv<-apply(DIFF,1,max)
LowerEnv<-apply(DIFF,1,min)
# Plot the data
plot(s, khcases-khcontrols, type="l", ylim=c(-.015,0.015))
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
plot(s, khcases-khcontrols, type="l", ylim=c(-.015,0.015))
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
plot(s, khcases-khcontrols, type="l")
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
UppEnv
sG<-Sobj_SpatialGrid(StudyArea)$SG
gt<-slot(sG,"grid")
# compute the risk ratio by estimating the intensity of cases and controls and taking the ratio
grid<-slot(slot(slot(StudyArea, "polygons")[[1]], "Polygons")[[1]], "coords")
mean(dist(coordinates(cases)))
s<-seq(0,0.12, by=0.0001) #A vector of distances at which to calculate the K function
# Calculate the K function for cases and controls over the grid
khcases<-khat(coordinates(cases), grid, s)
khcontrols<-khat(coordinates(controls), grid, s)
# Calculate the covariance matrix for the difference between two K-functions under random labelling of the corresponding two sets of points.
khcov<-khvmat(coordinates(cases), coordinates(controls), grid, s)
# Test statistic
T0<-sum(((khcases-khcontrols))/sqrt(diag(khcov+0.00000000001)))
niter<-999  # Number of iterations to simulate
T<-rep(NA,niter) # Empty vector in which to store simulated values of T
khcontrols
T0
niter<-999  # Number of iterations to simulate
T<-rep(NA,niter) # Empty vector in which to store simulated values of T
khcasesrel<-matrix(NA, nrow=length(s), ncol=niter)
khcontrolsrel<-matrix(NA, nrow=length(s), ncol=niter)
DIFF<-NULL
for(i in 1:niter){
idxrel<-sample(CaseControl$IndexCaseHousehold)==0
casesrel<-coordinates(CaseControl[idxrel,])
controlsrel<-coordinates(CaseControl[!idxrel,])
khcasesrel[,i]<-khat(casesrel, grid, s)
khcontrolsrel[,i]<-khat(controlsrel, grid, s)
khdiff<-khcasesrel[,i]-khcontrolsrel[,i]
DIFF<-cbind(DIFF,khdiff)
length(khdiff); length(sqrt(diag(khcov)))
T[i]<-sum(khdiff/sqrt(diag(khcov+0.00000000001)))  # check this line, did not work as entered
}
pvalue<-(sum(T>T0)+1)/(niter+1)
pvalue
khcases
khcontrols
khcasesrel
khcasesrel<-matrix(NA, nrow=length(s), ncol=niter)
khcontrolsrel<-matrix(NA, nrow=length(s), ncol=niter)
DIFF
DIFF<-NULL
i<-1
idxrel<-sample(CaseControl$IndexCaseHousehold)==0
casesrel<-coordinates(CaseControl[idxrel,])
controlsrel<-coordinates(CaseControl[!idxrel,])
controlsrel
idxrel
CaseControl$case
DIFF<-NULL
for(i in 1:niter){
idxrel<-sample(CaseControl$case)==0
casesrel<-coordinates(CaseControl[idxrel,])
controlsrel<-coordinates(CaseControl[!idxrel,])
khcasesrel[,i]<-khat(casesrel, grid, s)
khcontrolsrel[,i]<-khat(controlsrel, grid, s)
khdiff<-khcasesrel[,i]-khcontrolsrel[,i]
DIFF<-cbind(DIFF,khdiff)
length(khdiff); length(sqrt(diag(khcov)))
T[i]<-sum(khdiff/sqrt(diag(khcov+0.00000000001)))  # check this line, did not work as entered
}
# Calculate the P-value (proportion of times the simulated value is more extreme than observed)
pvalue<-(sum(T>T0)+1)/(niter+1)
pvalue
# Estimate 95% CI by taking the maximum and minimum simulated value for each point estimated
UppEnv<-apply(DIFF,1,max)
LowerEnv<-apply(DIFF,1,min)
# Plot the data
plot(s, khcases-khcontrols, type="l")
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
LowerEnv
UppEnv
plot(s, khcases-khcontrols, type="l", ylim=c(-1,1))
plot(s, khcases-khcontrols, type="l")
plot(s, khcases-khcontrols, type="l", ylim=c(-0.001,0.001))
plot(s, khcases-khcontrols, type="l", ylim=c(-0.001,0.01))
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
UppEnv<-apply(DIFF,1,quantile(x, .95) )
UppEnv<-apply(DIFF,1, quantile, probs = c(0.05, 0.9))
UppEnv<-apply(DIFF,1, quantile, probs = 0.95)
LowerEnv<-apply(DIFF,1,quantile, probs = 0.05)
# Plot the data
plot(s, khcases-khcontrols, type="l", ylim=c(-0.001,0.01))
lines(s, UppEnv, lty=2)
lines(s, LowerEnv, lty=2)
-s
polygon(c(s,-s),c(UppEnv, LowerEnv),col=gray(0.8))
polygon(c(s,s),c(UppEnv, LowerEnv),col=gray(0.8))
polygon(c(s,-s),c(UppEnv, -LowerEnv),col=gray(0.8))
polygon(c(s,rev(s)),c(UppEnv, rev(LowerEnv)),col=gray(0.8))
plot(s, khcases-khcontrols, type="l", ylim=c(-0.001,0.01))
polygon(c(s,rev(s)),c(UppEnv, rev(LowerEnv)),col=gray(0.8))
plot(s, khcases-khcontrols, type="l", ylim=c(-0.01,0.01))
polygon(c(s,rev(s)),c(UppEnv, rev(LowerEnv)),col=gray(0.8))
lines(s, khcases-khcontrols, type="l", col="red")
CaseControl<-read.csv(paste(root, path, data.dir,"CaseControl.csv", sep=""), header=T, stringsAsFactors=FALSE)
project2<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 " #
StudyArea<-readShapePoly("StudySiteDistricts.shp",IDvar=NULL, proj4string=CRS(project2))
# How many cases and controls are there?
sum(CaseControl$case, na.rm=T)
sum(CaseControl$case==0, na.rm=T)
# Export file for SatScan
cases<-CaseControl[which(CaseControl$case==1),]
controls<-CaseControl[which(CaseControl$case==0),]
cases$ncase<-1
controls$ncontrol<-1
cases<-cases[,c("household_id", "ncase")]
controls<-controls[,c("household_id", "ncontrol")]
locations<-CaseControl[,c("household_id", "lat", "long")]
